Mapping Areas That Could Benefit From Affordable Housing in New Jersey
Scott Mullarkey
04/29/2020
1.1 Abstract
Research has shown that affordable housing provides the community with a variety of benefits including improved housing affordability, neutral or positive impact on surrounding property values, an increase in diversity, reduced travel time to places of employment, an increase in municipal real estate tax collections, improved educational outcomes for children, lowered rates of asthma, and a reduction in poverty rates of seniors ages 65 and older. The State of New Jersey can use affordable housing development as a means of improving aspects of local communities by targeting its resources to areas that could benefit most from additional affordable housing.
A SHINY application was developed to highlight these areas. A user can select from a series of variables and the application with highlight in green which areas could benefit most from additional affordable housing based on the variables selected. The application draws from 10 years of data (2009-2018) for each variable and locations were scored based on their rate of change within that time period. The following R Markdown shows how this data was collected, processed and used within the application. Details on how the application was developed and can be used are also provided. The application can be accessed through the following link:
2.1 Introduction
This R Markdown was prepared to detail the development of a SHINY application that maps areas of New Jersey that could benefit most from additional affordable housing. It highlghts the data that was used as well as how it was processed using R. Code was provided for each segment. The document also highlights the application itself and how it is used. Code was provided for the development of the application.
3.1 Problem Statement
The intent of this project was to develop a SHINY application to highlight areas of New Jersey that could benefit most from additional affordable housing. This application would draw upon 10 years (2009-2018) of data on affordability, property values, diversity, travel time, real estate taxes, education, asthma, and poverty 65+ which would be processed in R. This application could be used by state officials to help determine the best means of resource allocation using affordable housing development as a tool to address a variety of social needs.
4.1 Literature Review
Research has shown that affordable housing provides the community with a variety of benefits. The following is a review of current literature that includes a brieft summary of each as weill as link to the original reports and publications.
4.1.1 Affordability
The National Low Income Housing Coalition publishes an annual report that analyzes the levels of housing affordability at the national, state, and county levels. Hourly wage needed to afford a two bedroom apartment - paying no more than 30% of their income towards housing - is one metric that is analyzed. Shortages of affordable rental housing leads to increases in this metric. (Source: National Low Income Housing Coalition)
4.1.2 Property Values
The Center for Housing Policy published a report that noted that affordable housing had a neutral or positive effect on the property values of neighboring parcels. (Source: Center for Housing Policy)
4.1.3 Diversity
The Brookings Institute and The Urban Institute published a report that noted production programs that increase availabilty of affordable housing promotes racial diversity. (Source: The Brookings Institute)
4.1.4 Travel Time
The Brookings Institute publised an article that noted an increase in affordable housing helps reduce commute times as provides lower-income workers increased housing options near places of employment. (Source: The Brookings Institute)
4.1.5 Real Estate Taxes
The Center for Housing Policy published a report that noted that affordable housing can lead to higher real estate income for a municipalit as it often reduces instances of foreclosure which leads to depressed property values and tax collections. (Source: Center for Housing Policy)
4.1.6 Education
The Center for Housing Policy published a report that noted affordable housing may improve educational outcomes of children as it helps increase housing stability reduces family homelessness. (Source: Center for Housing Policy)
4.1.7 Asthma
The American Public Health Assoication published an article that noted the green and sustainable development practices in affordable housing development and renovations can lead to lower asthma rates of residents. (Source: American Public Health Assoication)
4.1.8 Poverty 65+
The Urban Institute published a report that noted affordable housing helps reduce poverty in adults 65 and older as it allows them to age in place which is often a more affordable option than relocating to assisted-living facilities. (Source: The Urban Institute)
5.1 Methodology
The following methodology was used in collecting and processing the data as well as developing the application.
5.1.1 Application Design
An application was developed in SHINY to highlight areas of New Jersey that could benefit most from affordable housing. It draws upon 10 years (2009-2018) of data on affordability, property values, diversity, travel time, real estate taxes, education, asthma, and poverty 65+ which would be processed in R.
5.1.2 Data Preparation
A model was developed using R statistical package to aggregate and process all data to create a final csv file that will be used in the SHINY application. Details on data collection, analysis, processing, and aggregation are included below. Code for the SHINY application as well as a URL are also found below.
5.1.2.1 Setting Parameters
The model begins by downloading the required R libraries. It then sets up plot and map themes as well as a color palette. This allows the model to produce consistent visualizations throughout.
knitr::opts_chunk$set(cache=TRUE)
setwd("C:/Users/scott/Desktop/MUSA800_FINAL")
options(scipen=999)
library(tidyverse)
library(sf)
library(lubridate)
library(tigris)
library(tidycensus)
library(viridis)
library(riem)
library(gridExtra)
library(knitr)
library(kableExtra)
library(RSocrata)
library(gganimate)
library(ggplot2)
library(dplyr)
library(MASS)
library(QuantPsyc)
library(caret)
library(spatstat)
library(spdep)
library(FNN)
library(grid)
library(RColorBrewer)
library(wesanderson)
library(gridExtra)
library(grid)
library(stats)
library(lme4)
library(rgdal)
library(shiny)
plotTheme <- theme(
plot.title =element_text(size=12),
plot.subtitle = element_text(size=8),
plot.caption = element_text(size = 6),
axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
axis.text.y = element_text(size = 10),
axis.title.y = element_text(size = 10),
# Set the entire chart region to blank
panel.background=element_blank(),
plot.background=element_blank(),
#panel.border=element_rect(colour="#F0F0F0"),
# Format the grid
panel.grid.major=element_line(colour="#D0D0D0",size=.2),
axis.ticks=element_blank())
mapTheme <- function(base_size = 12) {
theme(
text = element_text( color = "black"),
plot.title = element_text(size = 14,colour = "black"),
plot.subtitle=element_text(face="italic"),
plot.caption=element_text(hjust=0),
axis.ticks = element_blank(),
panel.background = element_blank(),axis.title = element_blank(),
axis.text = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(colour = "black", fill=NA, size=2)
)
}
palette5 <- c("#eff3ff","#bdd7e7","#6baed6","#3182bd","#08519c")
palette4 <- c("#D2FBD4","#92BCAB","#527D82","#123F5A")
palette2 <- c("#6baed6","#08519c")5.1.2.2 Establishing Boundaries
State, County, and municipal boundaries are added to the model. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
#-------------------Create NJ Boundary-------------------
US_Boundary <- st_read("C:/Users/scott/Desktop/MUSA800/Data/states.shp")
US_Boundary <- US_Boundary %>% st_transform(crs=102271)
NJ_Boundary<- US_Boundary[(US_Boundary$STATE_ABBR=="NJ"),]
NJ_Boundary <-
NJ_Boundary %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 102271, agr = "constant")
NJ_Towns <- st_read("https://opendata.arcgis.com/datasets/3d5d1db8a1b34b418c331f4ce1fd0fef_2.geojson")
NJ_Towns <-
NJ_Towns %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 102271, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary))
NJ_Towns <- NJ_Towns %>% dplyr::select("NAME", "MUN_CODE", "geometry")
NJ_Counties <- st_read("https://opendata.arcgis.com/datasets/5f45e1ece6e14ef5866974a7b57d3b95_1.geojson")
NJ_Counties <-
NJ_Counties %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 102271, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary))
ggplot() +
geom_sf(data=NJ_Counties, fill=NA, colour="black") +
mapTheme()
5.1.2.3 Creating a Fishnet
A fishnet was applied to produce equal sized polygons (2,500 sqare meters) over the state and use them as a unit of measure for my analysis. Fishnet and boundaries were converted from crs of 102271 to 4326 so it would work in SHINY application.
knitr::opts_chunk$set(cache=TRUE)
#-------------------Create Fishnet-------------------
createFishnet <- function(x) {
fishnet <-
st_make_grid(NJ_Boundary, cellsize = x) %>%
st_sf()
fishnet <-
fishnet[NJ_Boundary,] %>%
mutate(uniqueID = rownames(.)) %>%
dplyr::select(uniqueID)
fishnet <-
fishnet %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 102271, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary))
fishnet <<- st_intersection(fishnet, NJ_Counties)%>%
dplyr::select(uniqueID, geometry,COUNTY)}
createFishnet(2500)
#-------------------Convert CRS to 4326 so it will work in SHINY-------------------
US_Boundary <- US_Boundary %>% st_transform(crs=4326)
NJ_Boundary<- US_Boundary[(US_Boundary$STATE_ABBR=="NJ"),]
NJ_Boundary <-
NJ_Boundary %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant")
NJ_Towns <- st_read("https://opendata.arcgis.com/datasets/3d5d1db8a1b34b418c331f4ce1fd0fef_2.geojson")
NJ_Towns <-
NJ_Towns %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary))
NJ_Towns <- NJ_Towns %>% dplyr::select("NAME", "MUN_CODE", "geometry")
NJ_Counties <- st_read("https://opendata.arcgis.com/datasets/5f45e1ece6e14ef5866974a7b57d3b95_1.geojson")
NJ_Counties <-
NJ_Counties %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary))
fishnet <-
fishnet %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary))
ggplot() +
geom_sf(data=NJ_Counties, fill=NA, colour="black") +
geom_sf(data=fishnet, fill=NA, colour="GRAY") +
labs(title = "Fishnet in New Jersey") +
mapTheme()
5.1.2.4 Adding Data
Ten years of data (2009-2018) was collected and processed for each variable as detailed below.
5.1.2.4.1 Data on Affordability
The model measures affordability by the hourly wage needed to afford a two bedroom apartment by county. Data was collected from the National Low Income Housing Coalition through their annual Out of Reach reports. County level data was pulled from reports of years 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 and entered into a csv file. This was then inputted into a dataframe with a seperate column for each year. The data was then merged with the fishnet so that every grid cell took on the underlying value in each year. Below is a facetted map series of the output.
knitr::opts_chunk$set(cache=TRUE)
#-------------------AFFORDABILITY-------------------
#Hourly wage necessary to afford 2 BR FMR
#https://nlihc.org/sites/default/files/oor/2009-OOR.pdf
#https://nlihc.org/sites/default/files/oor/2010-OOR.pdf
#https://nlihc.org/sites/default/files/oor/2011-OOR.pdf
#https://nlihc.org/sites/default/files/oor/2012-OOR.pdf
#https://reports.nlihc.org/sites/default/files/oor/2013-OOR-NJ_0.pdf
#https://reports.nlihc.org/sites/default/files/oor/2014-OOR-NJ.pdf
#https://nlihc.org/sites/default/files/oor/OOR_2015_FULL.pdf
#https://nlihc.org/sites/default/files/oor/OOR_2016.pdf
#https://reports.nlihc.org/sites/default/files/oor/OOR_2017_0.pdf
#https://nlihc.org/sites/default/files/oor/OOR_2018.pdf
NJ_Afford <- read.csv("C:/Users/scott/Desktop/MUSA800/Data/NJ_Affordability.csv")
names(NJ_Afford)[1] <- "COUNTY"
AFFORD.Net <- fishnet
AFFORD.Net$L2009 <- NJ_Afford$L2009[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)]
AFFORD.Net$L2010 <- NJ_Afford$L2010[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)]
AFFORD.Net$L2011 <- NJ_Afford$L2011[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)]
AFFORD.Net$L2012 <- NJ_Afford$L2012[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)]
AFFORD.Net$L2013 <- NJ_Afford$L2013[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)]
AFFORD.Net$L2014 <- NJ_Afford$L2014[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)]
AFFORD.Net$L2015 <- NJ_Afford$L2015[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)]
AFFORD.Net$L2016 <- NJ_Afford$L2016[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)]
AFFORD.Net$L2017 <- NJ_Afford$L2017[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)]
AFFORD.Net$L2018 <- NJ_Afford$L2018[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)]
##Facett Maps
AFFORD.Map <- AFFORD.Net %>%
gather(PIS, AFFORD, L2009:L2018)
AFFORD.Map$PIS <- substring(AFFORD.Map$PIS, 2)
ggplot() +
geom_sf(data = AFFORD.Map, aes(fill = AFFORD)) +
labs(title = "Wage Required-2 Bedroom Apartment, New Jersey") +
mapTheme() +
facet_wrap(~PIS, ncol = 5) +
scale_fill_gradient2("Value", low = "#762A83", mid = "white", high = "#1B7837")
The model then calculates the rate of change for each fishnet grid cell. This is done by looping an lm() function through the Affordability dataset which finds the slope of values for years 2009-2018 for each unique ID. The output shows the rate at which affordability has changed over time and which direction (increase/decrease) that change is moving. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
###Calculate Rate of Change
data <- AFFORD.Net %>%
dplyr::select(uniqueID,L2009:L2018)
data[is.na(data)] <- 0
data$geometry <- NULL
data$slope <- NA
data[1:nrow(data), "slope"] <- apply(data[1:nrow(data),
c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")],
1, function(row_i){lm(unlist(row_i) ~ unlist(data[1,
c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")]))$coefficients[2]
})
AFFORD.slope <- fishnet
AFFORD.slope$slope <- data$slope[match(data$uniqueID, AFFORD.slope$uniqueID)]
#RE.slope$slope[is.na(RE.slope$slope)] <- 0
ggplot() +
geom_sf(data = AFFORD.slope, aes(fill = slope)) +
labs(title = "Slope of Affordability (2009-2018), New Jersey") +
mapTheme() +
scale_fill_gradient2("Slope", low = "#762A83", mid = "white", high = "#1B7837")
A score is then applied to the slope values. This is done by assigning a value of “1” to the highest half of slope values which are those fishnet cells that have the fastest reduction in affordability, and a value of “2” for the lowest half of the slope values which are those fishnet cells that have the slowest reduction in affordability. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
###Assign Score
#https://stackoverflow.com/questions/7508229/how-to-create-a-column-with-a-quartile-rank/33896145
AFFORD.slope$slope2 <- AFFORD.slope$slope * -1
AFFORD.slope <- within(AFFORD.slope, score <- as.integer(cut(slope2, quantile(slope2, probs=0:2/2),
include.lowest=TRUE)))
ggplot() +
geom_sf(data = AFFORD.slope, aes(fill = score)) +
labs(title = "Affordabilty Score, New Jersey") +
mapTheme() +
scale_fill_gradient2("Score", limits = c(1, 2), labels=c("1","2"),breaks=c(1,2),low = "#762A83", high = "#1B7837")
5.1.2.4.2 Data on Property Value
The model measures property value by the mean property value by census tract. Data was collected from the US Census Bureau using tidycensus for a ten year period (2009-2018) saved to a dataframe with a seperate column for each year. The data was then merged with the fishnet so that every grid cell took on the average mean property value in each year. Here the average was taken to account for grid cells that had multiple census tract boudaries. Below is a facetted map series of the output.
knitr::opts_chunk$set(cache=TRUE)
#-------------------PROPERTY VALUE-------------------
#Enterprise Report: https://homeforallsmc.org/wp-content/uploads/2017/05/Impact-of-Affordable-Housing-on-Families-and-Communities.pdf
#Reference: https://furmancenter.org/files/media/Dont_Put_It_Here.pdf
#https://censusreporter.org/
census_api_key("7c4f7a4ee0b34622dcab1bef185348b87ae7b355", overwrite = TRUE)
NJCensus <-
get_acs(geography = "tract",
variables = c("B25077_001"),
year = 2010,
state = "NJ",
geometry = TRUE,
output = "wide") %>%
rename(PropVal = B25077_001E) %>%
dplyr::select(PropVal, GEOID, geometry)
NJCensus <-
NJCensus %>%
as.data.frame() %>%
distinct(GEOID, .keep_all = TRUE) %>%
#dplyr::select(GEOID, geometry) %>%
st_sf
NJCensus <-
NJCensus %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary),
join=st_intersects,
left = TRUE)
census <- function(x) {
NJCensus <<-
get_acs(geography = "tract",
variables = c("B25077_001"),
year = x,
state = "NJ",
geometry = TRUE,
output = "wide") %>%
rename(PropVal = B25077_001E) %>%
dplyr::select(PropVal, GEOID, geometry) %>%
as.data.frame() %>%
distinct(GEOID, .keep_all = TRUE) %>%
#dplyr::select(GEOID, geometry) %>%
st_sf%>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary),
join=st_intersects,
left = TRUE)
NJCensus$PIS <<-x
}
AVG <- function(x) {
NET <- NET[NET$PIS <= x,]
AVG <- aggregate( PropVal ~ uniqueID, NET, mean)
PROPVAL.Net$y <<- AVG$PropVal[match(PROPVAL.Net$uniqueID, AVG$uniqueID)]
PROPVAL.Net$y[is.na(PROPVAL.Net$y)] <<- 0
}
PROPVAL.Net <- fishnet
census(2010); NET <- st_intersection(fishnet, NJCensus)
AVG(2010); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2009"
census(2010); NET <- st_intersection(fishnet, NJCensus)
AVG(2010); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2010"
census(2011); NET <- st_intersection(fishnet, NJCensus)
AVG(2011); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2011"
census(2012); NET <- st_intersection(fishnet, NJCensus)
AVG(2012); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2012"
census(2012); NET <- st_intersection(fishnet, NJCensus)
AVG(2012); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2013"
census(2014); NET <- st_intersection(fishnet, NJCensus)
AVG(2014); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2014"
census(2015); NET <- st_intersection(fishnet, NJCensus)
AVG(2015); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2015"
census(2015); NET <- st_intersection(fishnet, NJCensus)
AVG(2015); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2016"
census(2017); NET <- st_intersection(fishnet, NJCensus)
AVG(2017); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2017"
census(2018); NET <- st_intersection(fishnet, NJCensus)
AVG(2018); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2018"
###Facett Maps
PROPVAL.Map <- PROPVAL.Net %>%
gather(PIS, PROPVAL, L2009:L2018)
PROPVAL.Map$PIS <- substring(PROPVAL.Map$PIS, 2)
ggplot() +
geom_sf(data = PROPVAL.Map, aes(fill = PROPVAL)) +
labs(title = "Property Value, New Jersey") +
mapTheme() +
facet_wrap(~PIS, ncol = 5) +
scale_fill_gradient2("Value", low = "#762A83", mid = "white", high = "#1B7837")
The model then calculates the rate of change for each fishnet grid cell. This is done by looping an lm() function through the Property Value dataset which finds the slope of values for years 2009-2018 for each unique ID. The output shows the rate at which mean property value has changed over time and which direction (increase/decrease) that change is moving. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
###Calculate Rate of Change
data <- PROPVAL.Net %>%
dplyr::select(uniqueID,L2009:L2018)
data$geometry <- NULL
data$slope <- NA
data[1:nrow(data), "slope"] <- apply(data[1:nrow(data),
c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")],
1, function(row_i){lm(unlist(row_i) ~ unlist(data[1,
c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")]))$coefficients[2]
})
PROPVAL.slope <- fishnet
PROPVAL.slope$slope <- data$slope[match(data$uniqueID, PROPVAL.slope$uniqueID)]
#PROPVAL.slope$slope[is.na(PROPVAL.slope$slope)] <- 0
ggplot() +
geom_sf(data = PROPVAL.slope, aes(fill = slope)) +
labs(title = "Slope of Property Value (2009-2018), New Jersey") +
mapTheme() +
scale_fill_gradient2("Slope", low = "#762A83", mid = "white", high = "#1B7837")
A score is then applied to the slope values. This is done by assigning a value of “1” to the lowest half of slope values which are those fishnet cells that have the slowest increase in property values, and a value of “2” for the highest half of the slope values which are those fishnet cells that have the fastest increase in property values. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
###Assign Score
#https://stackoverflow.com/questions/7508229/how-to-create-a-column-with-a-quartile-rank/33896145
#PROPVAL.slope$slope2 <- PROPVAL.slope$slope * -1
PROPVAL.slope <- within(PROPVAL.slope, score <- as.integer(cut(slope, quantile(slope, probs=0:2/2),
include.lowest=TRUE)))
ggplot() +
geom_sf(data = PROPVAL.slope, aes(fill = score)) +
labs(title = "Property Value Score, New Jersey") +
mapTheme() +
scale_fill_gradient2("Score", limits = c(1, 2), labels=c("1","2"),breaks=c(1,2),low = "#762A83", high = "#1B7837")
5.1.2.4.3 Data on Diversity
The model measures diversity by the percent white of population by census tract. Data was collected from the US Census Bureau using tidycensus for a ten year period (2009-2018) saved to a dataframe with a seperate column for each year. The data was then merged with the fishnet so that every grid cell took on the average percent white in each year. Here the average was taken to account for grid cells that had multiple census tract boudaries. Below is a facetted map series of the output.
knitr::opts_chunk$set(cache=TRUE)
#-------------------PERCENT WHITE-------------------
#Enterprise Report: https://homeforallsmc.org/wp-content/uploads/2017/05/Impact-of-Affordable-Housing-on-Families-and-Communities.pdf
#Reference: https://www.brookings.edu/wp-content/uploads/2016/06/housingreview.pdf
#https://censusreporter.org/
NJCensus <-
get_acs(geography = "tract",
variables = c("B01003_001","B02001_002"),
year = 2011,
state = "NJ",
geometry = TRUE,
output = "wide") %>%
rename(Total_Pop = B01003_001E, White_Pop = B02001_002E) %>%
dplyr::select(Total_Pop, White_Pop, GEOID, geometry) %>%
mutate(Percent_White = White_Pop / Total_Pop)
NJCensus <-
NJCensus %>%
as.data.frame() %>%
distinct(GEOID, .keep_all = TRUE) %>%
#dplyr::select(GEOID, geometry) %>%
st_sf
NJCensus <-
NJCensus %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary),
join=st_intersects,
left = TRUE)
census <- function(x) {
NJCensus <<-
get_acs(geography = "tract",
variables = c("B01003_001","B02001_002"),
year = x,
state = "NJ",
geometry = TRUE,
output = "wide") %>%
rename(Total_Pop = B01003_001E, White_Pop = B02001_002E) %>%
dplyr::select(Total_Pop, White_Pop, GEOID, geometry) %>%
mutate(Percent_White = White_Pop / Total_Pop)%>%
as.data.frame() %>%
distinct(GEOID, .keep_all = TRUE) %>%
#dplyr::select(GEOID, geometry) %>%
st_sf%>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary),
join=st_intersects,
left = TRUE)
NJCensus$PIS <<-x
}
AVG <- function(x) {
NET <- NET[NET$PIS <= x,]
AVG <- aggregate( Percent_White ~ uniqueID, NET, mean)
PCTWHITE.Net$y <<- AVG$Percent_White[match(PCTWHITE.Net$uniqueID, AVG$uniqueID)]
PCTWHITE.Net$y[is.na(PCTWHITE.Net$y)] <<- 0
}
PCTWHITE.Net <- fishnet
census(2010); NET <- st_intersection(fishnet, NJCensus)
AVG(2010); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2009"
census(2010); NET <- st_intersection(fishnet, NJCensus)
AVG(2010); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2010"
census(2011); NET <- st_intersection(fishnet, NJCensus)
AVG(2011); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2011"
census(2012); NET <- st_intersection(fishnet, NJCensus)
AVG(2012); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2012"
census(2012); NET <- st_intersection(fishnet, NJCensus)
AVG(2012); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2013"
census(2014); NET <- st_intersection(fishnet, NJCensus)
AVG(2014); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2014"
census(2015); NET <- st_intersection(fishnet, NJCensus)
AVG(2015); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2015"
census(2016); NET <- st_intersection(fishnet, NJCensus)
AVG(2016); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2016"
census(2017); NET <- st_intersection(fishnet, NJCensus)
AVG(2017); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2017"
census(2018); NET <- st_intersection(fishnet, NJCensus)
AVG(2018); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2018"
##Facett Maps
PCTWHITE.Map <- PCTWHITE.Net %>%
gather(PIS, PCTWHITE, L2009:L2018)
PCTWHITE.Map$PIS <- substring(PCTWHITE.Map$PIS, 2)
ggplot() +
geom_sf(data = PCTWHITE.Map, aes(fill = PCTWHITE)) +
labs(title = "Percent White, New Jersey") +
mapTheme() +
facet_wrap(~PIS, ncol = 5) +
scale_fill_gradient2("Value", low = "#762A83", mid = "white", high = "#1B7837")
The model then calculates the rate of change for each fishnet grid cell. This is done by looping an lm() function through the Diversity dataset which finds the slope of values for years 2009-2018 for each unique ID. The output shows the rate at which percent white has changed over time and which direction (increase/decrease) that change is moving. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
###Calculate Rate of Change
data <- PCTWHITE.Net %>%
dplyr::select(uniqueID,L2009:L2018)
data$geometry <- NULL
data$slope <- NA
data[1:nrow(data), "slope"] <- apply(data[1:nrow(data),
c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")],
1, function(row_i){lm(unlist(row_i) ~ unlist(data[1,
c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")]))$coefficients[2]
})
PCTWHITE.slope <- fishnet
PCTWHITE.slope$slope <- data$slope[match(data$uniqueID, PCTWHITE.slope$uniqueID)]
#PCTWHITE.slope$slope[is.na(PCTWHITE.slope$slope)] <- 0
ggplot() +
geom_sf(data = PCTWHITE.slope, aes(fill = slope)) +
labs(title = "Slope of Diversity (2009-2018), New Jersey") +
mapTheme() +
scale_fill_gradient2("Slope", low = "#762A83", mid = "white", high = "#1B7837")
A score is then applied to the slope values. This is done by assigning a value of “1” to the highest half of slope values which are those fishnet cells that have the fastest increase in percent white, and a value of “2” for the lowest half of the slope values which are those fishnet cells that have the fastest decrease in percent white. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
###Assign Score
#https://stackoverflow.com/questions/7508229/how-to-create-a-column-with-a-quartile-rank/33896145
PCTWHITE.slope$slope2 <- PCTWHITE.slope$slope * -1
PCTWHITE.slope <- within(PCTWHITE.slope, score <- as.integer(cut(slope2, quantile(slope2, probs=0:2/2),
include.lowest=TRUE)))
ggplot() +
geom_sf(data = PCTWHITE.slope, aes(fill = score)) +
labs(title = "Diversity Score, New Jersey") +
mapTheme() +
scale_fill_gradient2("Score", limits = c(1, 2), labels=c("1","2"),breaks=c(1,2),low = "#762A83", high = "#1B7837")
5.1.2.4.4 Data on Travel Time
The model measures travel time by percent of population that has over 90 minute commute to work by census tract. Data was collected from the US Census Bureau using tidycensus for a ten year period (2009-2018) saved to a dataframe with a seperate column for each year. The data was then merged with the fishnet so that every grid cell took on the average travel time in each year. Here the average was taken to account for grid cells that had multiple census tract boudaries. Below is a facetted map series of the output.
knitr::opts_chunk$set(cache=TRUE)
#-------------------TRAVEL TIME-------------------
#Enterprise Report: https://homeforallsmc.org/wp-content/uploads/2017/05/Impact-of-Affordable-Housing-on-Families-and-Communities.pdf
#Reference: http://www.reconnectingamerica.org/assets/Uploads/20120904AHpolicybrief.pdf
#https://censusreporter.org/
NJCensus <-
get_acs(geography = "tract",
variables = c("B01003_001","B08303_013"),
year = 2009,
state = "NJ",
geometry = TRUE,
output = "wide") %>%
rename(Total_Pop = B01003_001E, TRAVEL = B08303_013E) %>%
dplyr::select(Total_Pop, TRAVEL, GEOID, geometry) %>%
mutate(AVGCOMMUTE = TRAVEL / Total_Pop)
NJCensus <-
NJCensus %>%
as.data.frame() %>%
distinct(GEOID, .keep_all = TRUE) %>%
#dplyr::select(GEOID, geometry) %>%
st_sf
NJCensus <-
NJCensus %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary),
join=st_intersects,
left = TRUE)
census <- function(x) {
NJCensus <<-
get_acs(geography = "tract",
variables = c("B01003_001","B08303_013"),
year = x,
state = "NJ",
geometry = TRUE,
output = "wide") %>%
rename(Total_Pop = B01003_001E, TRAVEL = B08303_013E) %>%
dplyr::select(Total_Pop, TRAVEL, GEOID, geometry) %>%
mutate(AVGCOMMUTE = TRAVEL / Total_Pop) %>%
as.data.frame() %>%
distinct(GEOID, .keep_all = TRUE) %>%
#dplyr::select(GEOID, geometry) %>%
st_sf%>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary),
join=st_intersects,
left = TRUE)
NJCensus$PIS <<-x
}
AVG <- function(x) {
NET <- NET[NET$PIS <= x,]
AVG <- aggregate( AVGCOMMUTE ~ uniqueID, NET, mean)
TRAVEL.Net$y <<- AVG$AVGCOMMUTE[match(TRAVEL.Net$uniqueID, AVG$uniqueID)]
TRAVEL.Net$y[is.na(TRAVEL.Net$y)] <<- 0
}
TRAVEL.Net <- fishnet
census(2009); NET <- st_intersection(fishnet, NJCensus)
AVG(2009); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2009"
census(2010); NET <- st_intersection(fishnet, NJCensus)
AVG(2010); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2010"
census(2011); NET <- st_intersection(fishnet, NJCensus)
AVG(2011); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2011"
census(2012); NET <- st_intersection(fishnet, NJCensus)
AVG(2012); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2012"
census(2013); NET <- st_intersection(fishnet, NJCensus)
AVG(2013); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2013"
census(2014); NET <- st_intersection(fishnet, NJCensus)
AVG(2014); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2014"
census(2015); NET <- st_intersection(fishnet, NJCensus)
AVG(2015); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2015"
census(2015); NET <- st_intersection(fishnet, NJCensus)
AVG(2015); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2016"
census(2017); NET <- st_intersection(fishnet, NJCensus)
AVG(2017); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2017"
census(2018); NET <- st_intersection(fishnet, NJCensus)
AVG(2018); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2018"
##Facett Maps
TRAVEL.Map <- TRAVEL.Net %>%
gather(PIS, TRAVEL, L2009:L2018)
TRAVEL.Map$PIS <- substring(TRAVEL.Map$PIS, 2)
ggplot() +
geom_sf(data = TRAVEL.Map, aes(fill = TRAVEL)) +
labs(title = "Travel Time, New Jersey") +
mapTheme() +
facet_wrap(~PIS, ncol = 5) +
scale_fill_gradient2("Value", low = "#762A83", mid = "white", high = "#1B7837")
The model then calculates the rate of change for each fishnet grid cell. This is done by looping an lm() function through the Travel Time dataset which finds the slope of values for years 2009-2018 for each unique ID. The output shows the rate at which travel time has changed over time and which direction (increase/decrease) that change is moving. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
###Calculate Rate of Change
data <- TRAVEL.Net %>%
dplyr::select(uniqueID,L2009:L2018)
data$geometry <- NULL
data$slope <- NA
data[1:nrow(data), "slope"] <- apply(data[1:nrow(data),
c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")],
1, function(row_i){lm(unlist(row_i) ~ unlist(data[1,
c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")]))$coefficients[2]
})
TRAVEL.slope <- fishnet
TRAVEL.slope$slope <- data$slope[match(data$uniqueID, TRAVEL.slope$uniqueID)]
#TRAVEL.slope$slope[is.na(TRAVEL.slope$slope)] <- 0
ggplot() +
geom_sf(data = TRAVEL.slope, aes(fill = slope)) +
labs(title = "Slope of Travel Time (2009-2018), New Jersey") +
mapTheme() +
scale_fill_gradient2("Slope", low = "#762A83", mid = "white", high = "#1B7837")
A score is then applied to the slope values. This is done by assigning a value of “1” to the highest half of slope values which are those fishnet cells that have the fastest increase in travel, and a value of “2” for the lowest half of the slope values which are those fishnet cells that have the fastest decrease in travel. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
###Assign Score
#https://stackoverflow.com/questions/7508229/how-to-create-a-column-with-a-quartile-rank/33896145
#TRAVEL.slope$slope2 <- TRAVEL.slope$slope * -1
TRAVEL.slope <- within(TRAVEL.slope, score <- as.integer(cut(slope, quantile(slope, probs=0:2/2),
include.lowest=TRUE)))
ggplot() +
geom_sf(data = TRAVEL.slope, aes(fill = score)) +
labs(title = "Travel Score, New Jersey") +
mapTheme() +
scale_fill_gradient2("Score", limits = c(1, 2), labels=c("1","2"),breaks=c(1,2),low = "#762A83", high = "#1B7837")
5.1.2.4.5 Data on Taxes
The model measures real estate tax by the annual real estate tax collection by municipality. Data was collected from the NJ Department of Community Affairs for a ten year period (2009-2018) and then saved to a csv file. This was then inputted into a dataframe with a seperate column for each year. The data was then merged with the fishnet so that every grid cell took on the average mean real estate tax collections in each year. Here the average was taken to account for grid cells that had multiple municipality boudaries. Below is a facetted map series of the output.
knitr::opts_chunk$set(cache=TRUE)
#-------------------RE TAXES-------------------
#Enterprise Report: https://homeforallsmc.org/wp-content/uploads/2017/05/Impact-of-Affordable-Housing-on-Families-and-Communities.pdf
#Reference: https://providencehousing.org/wp-content/uploads/2014/03/Housing-and-Economic-Development-Report-2011.pdf
#https://censusreporter.org/
#Reference: https://providencehousing.org/wp-content/uploads/2014/03/Housing-and-Economic-Development-Report-2011.pdf
##Data Source: https://www.state.nj.us/dca/divisions/dlgs/resources/property_tax.html
NJ_RETaxes <- read.csv("C:/Users/scott/Desktop/MUSA800/Data/RE_Taxes.csv")
NJ_RETaxes <- NJ_RETaxes %>% rename(NAME = Town)
NJ_RETaxes <- merge(NJ_RETaxes,NJ_Towns,by="NAME")
NJ_RETaxes$X <- NULL
NJ_RETaxes$X.1 <- NULL
NJ_RETaxes[1] <- NULL
NJ_RETaxes <-
NJ_RETaxes %>%
st_as_sf(crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary))
NJ_RETaxes.Net <- st_intersection(fishnet, NJ_RETaxes)
RE.Net <- fishnet
AggTax <- aggregate( L2009 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2009 <- AggTax$L2009[match(RE.Net$uniqueID, AggTax$uniqueID)]
AggTax <- aggregate( L2010 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2010 <- AggTax$L2010[match(RE.Net$uniqueID, AggTax$uniqueID)]
AggTax <- aggregate( L2011 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2011 <- AggTax$L2011[match(RE.Net$uniqueID, AggTax$uniqueID)]
AggTax <- aggregate( L2012 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2012 <- AggTax$L2012[match(RE.Net$uniqueID, AggTax$uniqueID)]
AggTax <- aggregate( L2013 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2013 <- AggTax$L2013[match(RE.Net$uniqueID, AggTax$uniqueID)]
AggTax <- aggregate( L2014 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2014 <- AggTax$L2014[match(RE.Net$uniqueID, AggTax$uniqueID)]
AggTax <- aggregate( L2015 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2015 <- AggTax$L2015[match(RE.Net$uniqueID, AggTax$uniqueID)]
AggTax <- aggregate( L2016 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2016 <- AggTax$L2016[match(RE.Net$uniqueID, AggTax$uniqueID)]
AggTax <- aggregate( L2017 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2017 <- AggTax$L2017[match(RE.Net$uniqueID, AggTax$uniqueID)]
AggTax <- aggregate( L2018 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2018 <- AggTax$L2018[match(RE.Net$uniqueID, AggTax$uniqueID)]
##Facett Maps
RE.Map <- RE.Net %>%
gather(PIS, TAXES, L2009:L2018)
RE.Map$PIS <- substring(RE.Map$PIS, 2)
ggplot() +
geom_sf(data = RE.Map, aes(fill = TAXES)) +
labs(title = "Real Estate Taxes, New Jersey") +
mapTheme() +
facet_wrap(~PIS, ncol = 5) +
scale_fill_gradient2("Value", low = "#762A83", mid = "white", high = "#1B7837")
The model then calculates the rate of change for each fishnet grid cell. This is done by looping an lm() function through the Real Estate Tax dataset which finds the slope of values for years 2009-2018 for each unique ID. The output shows the rate at which real estate tax collection has changed over time and which direction (increase/decrease) that change is moving. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
###Calculate Rate of Change
data <- RE.Net %>%
dplyr::select(uniqueID,L2009:L2018)
data[is.na(data)] <- 0
data$geometry <- NULL
data$slope <- NA
data[1:nrow(data), "slope"] <- apply(data[1:nrow(data),
c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")],
1, function(row_i){lm(unlist(row_i) ~ unlist(data[1,
c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")]))$coefficients[2]
})
RE.slope <- fishnet
RE.slope$slope <- data$slope[match(data$uniqueID, RE.slope$uniqueID)]
#RE.slope$slope[is.na(RE.slope$slope)] <- 0
ggplot() +
geom_sf(data = RE.slope, aes(fill = slope)) +
labs(title = "Slope of Real Estate Taxes (2009-2018), New Jersey") +
mapTheme() +
scale_fill_gradient2("Slope", low = "#762A83", mid = "white", high = "#1B7837")
A score is then applied to the slope values. This is done by assigning a value of “1” to the lowest half of slope values which are those fishnet cells that have the slowest increase in real estate tax collections, and a value of “2” for the highest half of the slope values which are those fishnet cells that have the fastest increase in real estate tax collections. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
###Assign Score
#https://stackoverflow.com/questions/7508229/how-to-create-a-column-with-a-quartile-rank/33896145
#RE.slope$slope2 <- RE.slope$slope * -1
RE.slope <- within(RE.slope, score <- as.integer(cut(slope, quantile(slope, probs=0:2/2),
include.lowest=TRUE)))
ggplot() +
geom_sf(data = RE.slope, aes(fill = score)) +
labs(title = "Real Estate Tax Score, New Jersey") +
mapTheme() +
scale_fill_gradient2("Score", limits = c(1, 2), labels=c("1","2"),breaks=c(1,2),low = "#762A83", high = "#1B7837")
5.1.2.4.6 Data on Education
The model measures education by high school graduation by census tract. Data was collected from the US Census Bureau using tidycensus for a ten year period (2009-2018) saved to a dataframe with a seperate column for each year. The data was then merged with the fishnet so that every grid cell took on the average high school graduation in each year. Here the average was taken to account for grid cells that had multiple census tract boudaries. Below is a facetted map series of the output.
knitr::opts_chunk$set(cache=TRUE)
#-------------------EDUCATION-------------------
#Enterprise Report: https://homeforallsmc.org/wp-content/uploads/2017/05/Impact-of-Affordable-Housing-on-Families-and-Communities.pdf
#Reference: https://www.nhc.org/wp-content/uploads/2017/03/The-Impacts-of-Affordable-Housing-on-Education-1.pdf
#https://censusreporter.org/
NJCensus <-
get_acs(geography = "tract",
variables = c("B23006_009"),
year = 2009,
state = "NJ",
geometry = TRUE,
output = "wide") %>%
rename(HIGHSCHOOL = B23006_009E) %>%
dplyr::select(HIGHSCHOOL, GEOID, geometry)
NJCensus <-
NJCensus %>%
as.data.frame() %>%
distinct(GEOID, .keep_all = TRUE) %>%
#dplyr::select(GEOID, geometry) %>%
st_sf
NJCensus <-
NJCensus %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary),
join=st_intersects,
left = TRUE)
census <- function(x) {
NJCensus <<-
get_acs(geography = "tract",
variables = c("B23006_009"),
year = x,
state = "NJ",
geometry = TRUE,
output = "wide") %>%
rename(HIGHSCHOOL = B23006_009E) %>%
dplyr::select(HIGHSCHOOL, GEOID, geometry) %>%
as.data.frame() %>%
distinct(GEOID, .keep_all = TRUE) %>%
#dplyr::select(GEOID, geometry) %>%
st_sf%>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary),
join=st_intersects,
left = TRUE)
NJCensus$PIS <<-x
}
AVG <- function(x) {
NET <- NET[NET$PIS <= x,]
AVG <- aggregate( HIGHSCHOOL ~ uniqueID, NET, mean)
HIGHSCHOOL.Net$y <<- AVG$HIGHSCHOOL[match(HIGHSCHOOL.Net$uniqueID, AVG$uniqueID)]
HIGHSCHOOL.Net$y[is.na(HIGHSCHOOL.Net$y)] <<- 0
}
HIGHSCHOOL.Net <- fishnet
census(2009); NET <- st_intersection(fishnet, NJCensus)
AVG(2009); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2009"
census(2010); NET <- st_intersection(fishnet, NJCensus)
AVG(2010); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2010"
census(2011); NET <- st_intersection(fishnet, NJCensus)
AVG(2011); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2011"
census(2012); NET <- st_intersection(fishnet, NJCensus)
AVG(2012); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2012"
census(2013); NET <- st_intersection(fishnet, NJCensus)
AVG(2013); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2013"
census(2014); NET <- st_intersection(fishnet, NJCensus)
AVG(2014); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2014"
census(2015); NET <- st_intersection(fishnet, NJCensus)
AVG(2015); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2015"
census(2016); NET <- st_intersection(fishnet, NJCensus)
AVG(2016); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2016"
census(2017); NET <- st_intersection(fishnet, NJCensus)
AVG(2017); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2017"
census(2018); NET <- st_intersection(fishnet, NJCensus)
AVG(2018); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2018"
###Facett Maps
HIGHSCHOOL.Map <- HIGHSCHOOL.Net %>%
gather(PIS, HIGHSCHOOL, L2009:L2018)
HIGHSCHOOL.Map$PIS <- substring(HIGHSCHOOL.Map$PIS, 2)
ggplot() +
geom_sf(data = HIGHSCHOOL.Map, aes(fill = HIGHSCHOOL)) +
labs(title = "High School Degree, New Jersey") +
mapTheme() +
facet_wrap(~PIS, ncol = 5) +
scale_fill_gradient2("Value", low = "#762A83", mid = "white", high = "#1B7837")
The model then calculates the rate of change for each fishnet grid cell. This is done by looping an lm() function through the Education dataset which finds the slope of values for years 2009-2018 for each unique ID. The output shows the rate at which high school graduation has changed over time and which direction (increase/decrease) that change is moving. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
###Calculate Rate of Change
data <- HIGHSCHOOL.Net %>%
dplyr::select(uniqueID,L2009:L2018)
data$geometry <- NULL
data$slope <- NA
data[1:nrow(data), "slope"] <- apply(data[1:nrow(data),
c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")],
1, function(row_i){lm(unlist(row_i) ~ unlist(data[1,
c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")]))$coefficients[2]
})
HIGHSCHOOL.slope <- fishnet
HIGHSCHOOL.slope$slope <- data$slope[match(data$uniqueID, HIGHSCHOOL.slope$uniqueID)]
#PROPVAL.slope$slope[is.na(PROPVAL.slope$slope)] <- 0
ggplot() +
geom_sf(data = HIGHSCHOOL.slope, aes(fill = slope)) +
labs(title = "Slope of Education (2009-2018), New Jersey") +
mapTheme() +
scale_fill_gradient2("Slope", low = "#762A83", mid = "white", high = "#1B7837")
A score is then applied to the slope values. This is done by assigning a value of “1” to the lowest half of slope values which are those fishnet cells that have the slowest increase in high school graduation, and a value of “2” for the highest half of the slope values which are those fishnet cells that have the fastest increase in high school graduation. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
###Assign Score
#https://stackoverflow.com/questions/7508229/how-to-create-a-column-with-a-quartile-rank/33896145
#HIGHSCHOOL.slope$slope2 <- HIGHSCHOOL.slope$slope * -1
HIGHSCHOOL.slope <- within(HIGHSCHOOL.slope, score <- as.integer(cut(slope, quantile(slope, probs=0:2/2),
include.lowest=TRUE)))
ggplot() +
geom_sf(data = HIGHSCHOOL.slope, aes(fill = score)) +
labs(title = "Education Score, New Jersey") +
mapTheme() +
scale_fill_gradient2("Score", limits = c(1, 2), labels=c("1","2"),breaks=c(1,2),low = "#762A83", high = "#1B7837")
5.1.2.4.7 Data on Asthma
The model Asthma by Asthma hospitalizations by County. Data was collected from the NJ Department of Health for a ten year period (2009-2018) and saved to a csv file. This data was then inputted into a dataframe with a seperate column for each year. The data was then merged with the fishnet so that every grid cell took on the average travel time in each year. Here the average was taken to account for grid cells that had multiple census tract boudaries. Below is a facetted map series of the output.
knitr::opts_chunk$set(cache=TRUE)
#-------------------ASTHMA-------------------
#Enterprise Report: https://homeforallsmc.org/wp-content/uploads/2017/05/Impact-of-Affordable-Housing-on-Families-and-Communities.pdf
#Reference: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3000722/
#Source:https://www.nj.gov/health/fhs/chronic/documents/asthma_profiles/asthma_overview.pdf
#Source: https://www-doh.state.nj.us/doh-shad/query/selection/ub/UBSelection.html
#Age-adjusted Rates (Discharges Per 10,000 Standard Population) for Counties-AMI -- Primary Diagnosis
NJ_Asthma <- read.csv("C:/Users/scott/Desktop/MUSA800/Data/NJ_Asthma.csv")
names(NJ_Asthma)[1] <- "COUNTY"
ASTHMA.Net <- fishnet
ASTHMA.Net$L2009 <- NJ_Asthma$L2009[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)]
ASTHMA.Net$L2010 <- NJ_Asthma$L2010[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)]
ASTHMA.Net$L2011 <- NJ_Asthma$L2011[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)]
ASTHMA.Net$L2012 <- NJ_Asthma$L2012[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)]
ASTHMA.Net$L2013 <- NJ_Asthma$L2013[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)]
ASTHMA.Net$L2014 <- NJ_Asthma$L2014[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)]
ASTHMA.Net$L2015 <- NJ_Asthma$L2015[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)]
ASTHMA.Net$L2016 <- NJ_Asthma$L2016[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)]
ASTHMA.Net$L2017 <- NJ_Asthma$L2017[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)]
ASTHMA.Net$L2018 <- NJ_Asthma$L2018[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)]
##Facett Maps
ASTHMA.Map <- ASTHMA.Net %>%
gather(PIS, ASTHMA, L2009:L2018)
ASTHMA.Map$PIS <- substring(ASTHMA.Map$PIS, 2)
ggplot() +
geom_sf(data = ASTHMA.Map, aes(fill = ASTHMA)) +
labs(title = "Asthma Hospitalizations, New Jersey") +
mapTheme() +
facet_wrap(~PIS, ncol = 5) +
scale_fill_gradient2("Value", low = "#762A83", mid = "white", high = "#1B7837")
The model then calculates the rate of change for each fishnet grid cell. This is done by looping an lm() function through the Asthma dataset which finds the slope of values for years 2009-2018 for each unique ID. The output shows the rate at which asthma hospitalizations has changed over time and which direction (increase/decrease) that change is moving. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
###Calculate Rate of Change
data <- ASTHMA.Net %>%
dplyr::select(uniqueID,L2009:L2018)
data[is.na(data)] <- 0
data$geometry <- NULL
data$slope <- NA
data[1:nrow(data), "slope"] <- apply(data[1:nrow(data),
c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")],
1, function(row_i){lm(unlist(row_i) ~ unlist(data[1,
c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")]))$coefficients[2]
})
ASTHMA.slope <- fishnet
ASTHMA.slope$slope <- data$slope[match(data$uniqueID, RE.slope$uniqueID)]
#RE.slope$slope[is.na(RE.slope$slope)] <- 0
ggplot() +
geom_sf(data = ASTHMA.slope, aes(fill = slope)) +
labs(title = "Slope of Asthma Hospitalizations (2009-2018), New Jersey") +
mapTheme() +
scale_fill_gradient2("Slope", low = "#762A83", mid = "white", high = "#1B7837")
A score is then applied to the slope values. This is done by assigning a value of “1” to the highest half of slope values which are those fishnet cells that have the fastest increase in asthma hospitalizations, and a value of “2” for the lowest half of the slope values which are those fishnet cells that have the fastest decrease in asthma hospitalizations. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
###Assign Score
#https://stackoverflow.com/questions/7508229/how-to-create-a-column-with-a-quartile-rank/33896145
#ASTHMA.slope$slope2 <- ASTHMA.slope$slope * -1
ASTHMA.slope <- within(ASTHMA.slope, score <- as.integer(cut(slope, quantile(slope, probs=0:2/2),
include.lowest=TRUE)))
ggplot() +
geom_sf(data = ASTHMA.slope, aes(fill = score)) +
labs(title = "Asthma Score, New Jersey") +
mapTheme() +
scale_fill_gradient2("Score", limits = c(1, 2), labels=c("1","2"),breaks=c(1,2),low = "#762A83", high = "#1B7837")
5.1.2.4.8 Data on Poverty 65+
The model measures poverty 65+ by poverty status of households 65 and older by census tract. Data was collected from the US Census Bureau using tidycensus for a ten year period (2009-2018) saved to a dataframe with a seperate column for each year. The data was then merged with the fishnet so that every grid cell took on the poverty 65+ in each year. Here the average was taken to account for grid cells that had multiple census tract boudaries. Below is a facetted map series of the output.
knitr::opts_chunk$set(cache=TRUE)
#-------------------SENIORS-------------------
#Enterprise Report: https://homeforallsmc.org/wp-content/uploads/2017/05/Impact-of-Affordable-Housing-on-Families-and-Communities.pdf
#Reference: https://www.urban.org/research/publication/housing-platform-improving-outcomes-older-renters
#https://censusreporter.org/
NJCensus <-
get_acs(geography = "tract",
variables = c("B17017_008"),
year = 2009,
state = "NJ",
geometry = TRUE,
output = "wide") %>%
rename(POVERTY65 = B17017_008E) %>%
dplyr::select(POVERTY65, GEOID, geometry)
NJCensus <-
NJCensus %>%
as.data.frame() %>%
distinct(GEOID, .keep_all = TRUE) %>%
#dplyr::select(GEOID, geometry) %>%
st_sf
NJCensus <-
NJCensus %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary),
join=st_intersects,
left = TRUE)
census <- function(x) {
NJCensus <<-
get_acs(geography = "tract",
variables = c("B17017_008"),
year = x,
state = "NJ",
geometry = TRUE,
output = "wide") %>%
rename(POVERTY65 = B17017_008E) %>%
dplyr::select(POVERTY65, GEOID, geometry) %>%
as.data.frame() %>%
distinct(GEOID, .keep_all = TRUE) %>%
#dplyr::select(GEOID, geometry) %>%
st_sf%>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary),
join=st_intersects,
left = TRUE)
NJCensus$PIS <<-x
}
AVG <- function(x) {
NET <- NET[NET$PIS <= x,]
AVG <- aggregate( POVERTY65 ~ uniqueID, NET, mean)
POVERTY65.Net$y <<- AVG$POVERTY65[match(POVERTY65.Net$uniqueID, AVG$uniqueID)]
POVERTY65.Net$y[is.na(POVERTY65.Net$y)] <<- 0
}
POVERTY65.Net <- fishnet
census(2009); NET <- st_intersection(fishnet, NJCensus)
AVG(2009); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2009"
census(2010); NET <- st_intersection(fishnet, NJCensus)
AVG(2010); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2010"
census(2011); NET <- st_intersection(fishnet, NJCensus)
AVG(2011); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2011"
census(2012); NET <- st_intersection(fishnet, NJCensus)
AVG(2012); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2012"
census(2013); NET <- st_intersection(fishnet, NJCensus)
AVG(2013); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2013"
census(2014); NET <- st_intersection(fishnet, NJCensus)
AVG(2014); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2014"
census(2015); NET <- st_intersection(fishnet, NJCensus)
AVG(2015); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2015"
census(2016); NET <- st_intersection(fishnet, NJCensus)
AVG(2016); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2016"
census(2017); NET <- st_intersection(fishnet, NJCensus)
AVG(2017); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2017"
census(2018); NET <- st_intersection(fishnet, NJCensus)
AVG(2018); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2018"
###Facett Maps
POVERTY65.Map <- POVERTY65.Net %>%
gather(PIS, POVERTY65, L2009:L2018)
POVERTY65.Map$PIS <- substring(POVERTY65.Map$PIS, 2)
ggplot() +
geom_sf(data = POVERTY65.Map, aes(fill = POVERTY65)) +
labs(title = "Poverty Households 65+, New Jersey") +
mapTheme() +
facet_wrap(~PIS, ncol = 5) +
scale_fill_gradient2("Value", low = "#762A83", mid = "white", high = "#1B7837")
The model then calculates the rate of change for each fishnet grid cell. This is done by looping an lm() function through the Poverty 65+ dataset which finds the slope of values for years 2009-2018 for each unique ID. The output shows the rate at which poverty status 65+ has changed over time and which direction (increase/decrease) that change is moving. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
###Calculate Rate of Change
data <- POVERTY65.Net %>%
dplyr::select(uniqueID,L2009:L2018)
data$geometry <- NULL
data$slope <- NA
data[1:nrow(data), "slope"] <- apply(data[1:nrow(data),
c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")],
1, function(row_i){lm(unlist(row_i) ~ unlist(data[1,
c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")]))$coefficients[2]
})
POVERTY65.slope <- fishnet
POVERTY65.slope$slope <- data$slope[match(data$uniqueID, POVERTY65.slope$uniqueID)]
#PROPVAL.slope$slope[is.na(PROPVAL.slope$slope)] <- 0
ggplot() +
geom_sf(data = POVERTY65.slope, aes(fill = slope)) +
labs(title = "Slope of Poverty Households 65+ (2009-2018), New Jersey") +
mapTheme() +
scale_fill_gradient2("Slope", low = "#762A83", mid = "white", high = "#1B7837")
A score is then applied to the slope values. This is done by assigning a value of “1” to the highest half of slope values which are those fishnet cells that have the fastest increase in poverty status 65+, and a value of “2” for the lowest half of the slope values which are those fishnet cells that have the fastest decrease in poverty status 65+. Below is a map of the output.
knitr::opts_chunk$set(cache=TRUE)
###Assign Score
#https://stackoverflow.com/questions/7508229/how-to-create-a-column-with-a-quartile-rank/33896145
#POVERTY65.slope$slope2 <- POVERTY65.slope$slope * -1
POVERTY65.slope <- within(POVERTY65.slope, score <- as.integer(cut(slope, quantile(slope, probs=0:2/2),
include.lowest=TRUE)))
ggplot() +
geom_sf(data = POVERTY65.slope, aes(fill = score)) +
labs(title = "Poverty 65+ Score, New Jersey") +
mapTheme() +
scale_fill_gradient2("Score", limits = c(1, 2), labels=c("1","2"),breaks=c(1,2),low = "#762A83", high = "#1B7837")
5.1.3 Data Interpretation
All scores are combined into one dataframe. The dataframe is then set to all values greater than “1” to “0”. This leaves a file where all grid cells that could benefit from affordable housing have a value of “1” and all others have a value of “0”. The data is then saved to a csv file for use in the SHINY application. Below is a facetted map series of the output.
knitr::opts_chunk$set(cache=TRUE)
#Combine all scores to one dataframe
study.panel <- fishnet
study.panel$AFFORDABILITY <- AFFORD.slope$score[match(study.panel$uniqueID, PROPVAL.slope$uniqueID)]
study.panel$PROPERTYVALUE <- PROPVAL.slope$score[match(study.panel$uniqueID, PROPVAL.slope$uniqueID)]
study.panel$DIVERSITY <- PCTWHITE.slope$score[match(study.panel$uniqueID, PCTWHITE.slope$uniqueID)]
study.panel$TRAVELTIME <- TRAVEL.slope$score[match(study.panel$uniqueID, TRAVEL.slope$uniqueID)]
study.panel$RETAXES <- RE.slope$score[match(study.panel$uniqueID, RE.slope$uniqueID)]
study.panel$ASTHMA <- ASTHMA.slope$score[match(study.panel$uniqueID, ASTHMA.slope$uniqueID)]
study.panel$EDUCATION <- HIGHSCHOOL.slope$score[match(study.panel$uniqueID, HIGHSCHOOL.slope$uniqueID)]
study.panel$POVERTY65 <- POVERTY65.slope$score[match(study.panel$uniqueID, POVERTY65.slope$uniqueID)]
#Keep only scores of 1
study.panel$AFFORDABILITY[study.panel$AFFORDABILITY > 1] <- 0
study.panel$PROPERTYVALUE[study.panel$PROPERTYVALUE > 1] <- 0
study.panel$DIVERSITY[study.panel$DIVERSITY > 1] <- 0
study.panel$TRAVELTIME[study.panel$TRAVELTIME > 1] <- 0
study.panel$RETAXES[study.panel$RETAXES > 1] <- 0
study.panel$ASTHMA[study.panel$ASTHMA > 1] <- 0
study.panel$EDUCATION[study.panel$EDUCATION > 1] <- 0
study.panel$POVERTY65[study.panel$POVERTY65 > 1] <- 0
#Save as CSV and export to SHINY Application Folder
study.panel.csv <- study.panel
study.panel.csv$geometry <- NULL
write.csv(study.panel.csv,"C:/Users/scott/Desktop/MUSA800/SHINY/appData.csv", row.names = FALSE)
#Facet Map
study.panel.Map <- study.panel %>%
gather(Variable, Score, 4:11)
study.panel.Map$Score[study.panel.Map$Score == 0] <- NA
ggplot() +
geom_sf(data = study.panel.Map, aes(fill = Score)) +
labs(title = "Benefit from affordable housing") +
mapTheme() +
facet_wrap(~Variable, ncol = 4) +
scale_fill_gradientn(colours = "darkgreen")
#Export Fishnet to SHINY Application Folder
#setwd("C:/Users/scott/Desktop/MUSA800/SHINY")
#st_write(fishnet, "fishnet.shp")
#setwd("C:/Users/scott/Desktop/MUSA800")5.1.4 Data Presentation
An application was developed using SHINY. This data runs off the csv file that was created with each of the variable scores. Additional map details such as sample municipalities and major roads were included for reference.
knitr::opts_chunk$set(cache=TRUE)
#setwd("C:/Users/scott/Desktop/MUSA800/SHINY")
options(scipen=999)
library(shiny)
library(tidyverse)
library(sf)
library(lubridate)
library(tigris)
library(tidycensus)
library(viridis)
library(riem)
library(gridExtra)
library(knitr)
library(kableExtra)
library(RSocrata)
library(gganimate)
library(ggplot2)
library(dplyr)
library(MASS)
library(QuantPsyc)
library(caret)
library(spatstat)
library(spdep)
library(FNN)
library(grid)
library(RColorBrewer)
library(wesanderson)
library(gridExtra)
library(grid)
library(stats)
library(lme4)
library(rgdal)
library(sp)
library(ggrepel)
library(ggspatial)
plotTheme <- theme(
plot.title =element_text(size=12),
plot.subtitle = element_text(size=8),
plot.caption = element_text(size = 6),
axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
axis.text.y = element_text(size = 10),
axis.title.y = element_text(size = 10),
# Set the entire chart region to blank
panel.background=element_blank(),
plot.background=element_blank(),
#panel.border=element_rect(colour="#F0F0F0"),
# Format the grid
panel.grid.major=element_line(colour="#D0D0D0",size=.2),
axis.ticks=element_blank())
mapTheme <- function(base_size = 200) {
theme(
text = element_text( color = "black"),
plot.title = element_text(size = 14,colour = "black"),
plot.subtitle=element_text(face="italic"),
plot.caption=element_text(hjust=0),
axis.ticks = element_blank(),
panel.background = element_blank(),axis.title = element_blank(),
axis.text = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(colour = "black", fill=NA, size=2)
)
}
#-------------------Create NJ Boundary-------------------
US_Boundary <- st_read("states.shp")
US_Boundary <- US_Boundary %>% st_transform(crs=4326)
NJ_Boundary<- US_Boundary[(US_Boundary$STATE_ABBR=="NJ"),]
NJ_Boundary <-
NJ_Boundary %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant")
NJ_Counties <- st_read("https://opendata.arcgis.com/datasets/5f45e1ece6e14ef5866974a7b57d3b95_1.geojson")
NJ_Counties <-
NJ_Counties %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary))
#NJ_Counties = cbind(NJ_Counties, st_coordinates(st_centroid(NJ_Counties$geometry)))
NJ_Cities <- st_read("GISPORTAL_GISOWNER01_USCITIESTOWNS1MIL14.shp")
NJ_Cities <- st_intersection(NJ_Cities, NJ_Boundary)
NJ_Cities <- NJ_Cities %>% filter((NAME == "Mays Landing") | (NAME == "Hackensack") |
(NAME == "Mount Holly") | (NAME == "Camden") |
(NAME == "Cape May Court House") | (NAME == "Bridgeton") |
(NAME == "Newark") | (NAME == "Woodbury") |
(NAME == "Jersey City") | (NAME == "Flemington") |
(NAME == "New Brunswick") | (NAME == "Freehold") |
(NAME == "Morristown") | (NAME == "Toms River") |
(NAME == "Paterson") | (NAME == "Salem") |
(NAME == "Somerville") | (NAME == "Newton") |
(NAME == "Elizabeth") | (NAME == "Belvidere") )%>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary))
NJ_Roads <- st_read("tl_2015_34_prisecroads.shp")
GSP <- NJ_Roads %>% filter((FULLNAME == "Garden State Pkwy") )%>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary)) %>%
group_by(FULLNAME) %>% summarize()
GSP = cbind(GSP, st_coordinates(st_centroid(GSP$geometry)))
GSP$NAME <- "Garden State Parkway"
NJTPK <- NJ_Roads %>% filter((FULLNAME == "New Jersey Tpke") )%>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary)) %>%
group_by(FULLNAME) %>% summarize()
NJTPK = cbind(NJTPK, st_coordinates(st_centroid(NJTPK$geometry)))
NJTPK$NAME <- "NJ Turnpike"
I80 <- NJ_Roads %>% filter((FULLNAME == "I- 80") )%>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary)) %>%
group_by(FULLNAME) %>% summarize()
I80 = cbind(I80, st_coordinates(st_centroid(I80$geometry)))
I80$NAME <- "I-80"
I195 <- NJ_Roads %>% filter((FULLNAME == "I- 195") )%>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary)) %>%
group_by(FULLNAME) %>% summarize()
I195 = cbind(I195, st_coordinates(st_centroid(I195$geometry)))
I195$NAME <- "I-195"
US322 <- NJ_Roads %>% filter((FULLNAME == "US Hwy 322") )%>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary)) %>%
group_by(FULLNAME) %>% summarize()
US322 = cbind(US322, st_coordinates(st_centroid(US322$geometry)))
US322$NAME <- "US-322"
I78 <- NJ_Roads %>% filter((FULLNAME == "I- 78") )%>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary)) %>%
group_by(FULLNAME) %>% summarize()
I78 = cbind(I78, st_coordinates(st_centroid(I78$geometry)))
I78$NAME <- "I-78"
US206 <- NJ_Roads %>% filter((FULLNAME == "US Hwy 206") )%>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
st_transform(st_crs(NJ_Boundary)) %>%
group_by(FULLNAME) %>% summarize()
US206 = cbind(US206, st_coordinates(st_centroid(US206$geometry)))
US206$NAME <- "US-206"
#-------------------Create Fishnet-------------------
fishnet <- st_read("fishnet.shp")
fishnet <- st_intersection(fishnet, NJ_Counties)%>%
dplyr::select(uniqueID, geometry,COUNTY)
#-------------------Add Data-------------------
data <- read.csv("appData.csv")
study.panel <- fishnet
study.panel$AFFORDABILITY <- 1
study.panel$PROPERTYVALUE <- 1
study.panel$DIVERSITY <- 1
study.panel$TRAVELTIME <- 1
study.panel$RETAXES <- 1
study.panel$ASTHMA <- 1
study.panel$EDUCATION <- 1
study.panel$POVERTY65 <- 1
study.panel$TOTAL <- 1
#-------------------Create App-------------------
shinyApp(
ui = shinyUI(fluidPage(
titlePanel("Mapping Areas That Could Benefit Most From Affordable Housing in New Jersey"),
sidebarLayout(
sidebarPanel(
checkboxGroupInput("fun", label = "Variables",
choices = list("AFFORDABILITY" = 1,
"PROPERTY VALUE" = 2,
"DIVERSITY" = 3,
"TRAVEL TIME" = 4,
"REAL ESTATE TAXES" = 5,
"ASTHMA" = 6,
"EDUCATION" = 7,
"POVERTY 65+" = 8)
),
h5("Research has shown that affordable housing provides the community with a variety of benefits including improved housing affordability, neutral or positive impact on surrounding property values, an increase in diversity, reduced travel time to places of employment, an increase in municipal real estate tax collections, improved educational outcomes for children, lowered rates of asthma, and a reduction in poverty rates of seniors ages 65 and older."),
h5("This SHINY application was developed to highlight these areas. A user can select from the above variables and the application with highlight in green which areas could benefit most from additional affordable housing based on the variables selected."),
h5("The application draws from 10 years of data (2009-2018) for each variable and locations were scored based on their rate of change within that time period. The following R Markdown shows how this data was collected, processed and used within the application. Details on how the application was developed are also provided."),
url <- a("R Markdown", href="https://smullarkupenn.github.io/MUSA800_Final/ScottMullarkey_Final.html")
),
mainPanel(
plotOutput("plot")
)
)
)),
server = shinyServer(function(input, output) {
output$plot <- renderPlot({
if(1 %in% input$fun)
study.panel$AFFORDABILITY <- data$AFFORDABILITY[match(study.panel$uniqueID, data$uniqueID)]
if(2 %in% input$fun)
study.panel$PROPERTYVALUE <- data$PROPERTYVALUE[match(study.panel$uniqueID, data$uniqueID)]
if(3 %in% input$fun)
study.panel$DIVERSITY <- data$DIVERSITY[match(study.panel$uniqueID, data$uniqueID)]
if(4 %in% input$fun)
study.panel$TRAVELTIME <- data$TRAVELTIME[match(study.panel$uniqueID, data$uniqueID)]
if(5 %in% input$fun)
study.panel$RETAXES <- data$RETAXES[match(study.panel$uniqueID, data$uniqueID)]
if(6 %in% input$fun)
study.panel$ASTHMA <- data$ASTHMA[match(study.panel$uniqueID, data$uniqueID)]
if(7 %in% input$fun)
study.panel$EDUCATION <- data$EDUCATION[match(study.panel$uniqueID, data$uniqueID)]
if(8 %in% input$fun)
study.panel$POVERTY65 <- data$POVERTY65[match(study.panel$uniqueID, data$uniqueID)]
study.panel$TOTAL <- study.panel$AFFORDABILITY * study.panel$PROPERTYVALUE *
study.panel$DIVERSITY * study.panel$TRAVELTIME *
study.panel$RETAXES * study.panel$ASTHMA *
study.panel$EDUCATION * study.panel$POVERTY65
if(sum(study.panel$TOTAL)/ nrow(study.panel) == 1){study.panel$TOTAL <- 0}
study.panel$TOTAL[study.panel$TOTAL == 0] <- NA
ggplot() +
geom_sf(data = study.panel, aes(fill = TOTAL), color="WHITE") +
geom_sf(data=NJ_Counties, fill=NA, colour="BLACK") +
geom_sf(data=NJ_Cities, fill=NA, colour="BLUE") +
scale_fill_gradientn(colours = "GREEN", na.value="WHITE", guide=FALSE)+
geom_sf(data=GSP, fill=NA, colour="GOLD2") + geom_label(data=GSP,aes(x=X,y=Y,label=NAME), size=4) +
geom_sf(data=I80, fill=NA, colour="GOLD2") + geom_label(data=I80,aes(x=X,y=Y,label=NAME), size=4) +
geom_sf(data=I195, fill=NA, colour="GOLD2") + geom_label(data=I195,aes(x=X,y=Y,label=NAME), size=4) +
geom_sf(data=US322, fill=NA, colour="GOLD2") + geom_label(data=US322,aes(x=X,y=Y,label=NAME), size=4) +
geom_sf(data=I78, fill=NA, colour="GOLD2") + geom_label(data=I78,aes(x=X,y=Y,label=NAME), size=4) +
geom_sf(data=US206, fill=NA, colour="GOLD2") + geom_label(data=US206,aes(x=X,y=Y,label=NAME), size=4) +
geom_sf(data=NJTPK, fill=NA, colour="GOLD2") + geom_label(data=NJTPK,aes(x=X,y=Y,label=NAME), size=4) +
geom_text_repel(data = NJ_Cities, aes(LONGITUDE, LATITUDE, label = NAME), colour = "BLUE",
nudge_y = 0.05, size=3) +
annotation_north_arrow(location = "br", which_north = "true", style = north_arrow_fancy_orienteering)+
annotation_scale(location = "bl") +
mapTheme()
}, height = 750, width = 750)
})
)5.1.4.1 Accessing the Application
The application can be accessed through the following link: https://scottmullarkey.shinyapps.io/MUSA800_Final/
5.1.4.2 Demonstrating the Application
The application opens to the following page which shows a map of New Jersey with county boundaries, a variety of municiple locations for reference, as well as major highways. A description of the application is in the side bar on the left. A list of variables is included with a check box to their left
If the user selects one variable the map will highlight the areas of New Jersey that could benefit from additional affordable housing based on that one variable. Example is shown below. 
If the user selects multiple variables the map will highlight the areas of New Jersey that could benefit from additional affordable housing based on those multiple variables. Example is shown below. 
If the user selects all variables the map will highlight the areas of New Jersey that could benefit from additional affordable housing based on all variables. Example is shown below. 
6.1 Conclusion
This SHINY application was developed to highlight the areas that could benefit most from additional affordable housing based on current research. It can be used by the State of New Jersey to use affordable housing development as a means of improving aspects of local communities. A variety of users within an administration may be interested in such an application as it utilizes eight variables that cover a wide range of public interests.